# create subsets of data for manual RDD
smol1 <- data61 %>% 
  filter(lamp > -31 & lamp < 31)

smol2 <- data62 %>% 
  filter(lamp > -31 & lamp < 31)

smol3 <- data71 %>% 
  filter(one > -31 & one < 31)

smol4 <- data72 %>% 
  filter(one > -31 & one < 31)

smol5 <- data71 %>% 
  filter(two > -31 & two < 31)

smol6 <- data72 %>% 
  filter(two > -31 & two < 31)

smol7 <- data71 %>% 
  filter(three > -31 & three < 31)

smol8 <- data72 %>% 
  filter(three > -31 & three < 31)

smol9 <- data71 %>% 
  filter(four > -31 & four < 31)

smol10 <- data72 %>% 
  filter(four > -31 & four < 31)

smol11 <- data81 %>% 
  filter(five > -31 & five < 31)

smol12 <- data82 %>% 
  filter(five > -31 & five < 31)

smol13 <- data81 %>% 
  filter(six > -31 & six < 31)

smol14 <- data82 %>% 
  filter(six > -31 & six < 31)

smol15 <- data81 %>% 
  filter(seven > -31 & seven < 31)

smol16 <- data82 %>% 
  filter(seven > -31 & seven < 31)

smol17 <- data81 %>% 
  filter(eight > -31 & eight < 31)

smol18 <- data82 %>% 
  filter(eight > -31 & eight < 31)

smol19 <- data81 %>% 
  filter(nine > -31 & nine < 31)

smol20 <- data82 %>% 
  filter(nine > -31 & nine < 31)

smol21 <- data81 %>% 
  filter(ten > -31 & ten < 31)

smol22 <- data82 %>% 
  filter(ten > -31 & ten < 31)

smol23 <- data81 %>% 
  filter(eleven > -31 & eleven < 31)

smol24 <- data82 %>% 
  filter(eleven > -31 & eleven < 31)

smol25 <- data81 %>% 
  filter(twelve > -31 & twelve < 31)

smol26 <- data82 %>% 
  filter(twelve > -31 & twelve < 31)

smol27 <- data91 %>% 
  filter(thirteen > -31 & thirteen < 31)

smol28 <- data92 %>% 
  filter(thirteen > -31 & thirteen < 31)

smol29 <- data91IT %>% 
  filter(thirteen > -31 & thirteen < 31)

smol30 <- data92IT %>% 
  filter(thirteen > -31 & thirteen < 31)

smol31 <- data91 %>% 
  filter(fourteen > -31 & fourteen < 31)

smol32 <- data92 %>% 
  filter(fourteen > -31 & fourteen < 31)

smol33 <- data101 %>% 
  filter(fifteen > -31 & fifteen < 31)

smol34 <- data102 %>% 
  filter(fifteen > -31 & fifteen < 31)

smol35 <- data101IT %>% 
  filter(fifteen > -31 & fifteen < 31)

smol36 <- data102IT %>% 
  filter(fifteen > -31 & fifteen < 31)

smol37 <- CCnt %>% 
  filter(run > -31 & run < 31)

smol38 <- ITnt %>% 
  filter(run > -31 & run < 31)

smol39 <- CCon %>% 
  filter(run > -31 & run < 31)

smol40 <- ITon %>% 
  filter(run > -31 & run < 31)


# add treatment variable
smol1$tr <- 0
smol1$tr[smol1$lamp > -1] <- 1

smol2$tr <- 0
smol2$tr[smol2$lamp > -1] <- 1

smol3$tr <- 0
smol3$tr[smol3$one > -1] <- 1

smol4$tr <- 0
smol4$tr[smol4$one > -1] <- 1

smol5$tr <- 0
smol5$tr[smol5$two > -1] <- 1

smol6$tr <- 0
smol6$tr[smol4$two > -1] <- 1

smol7$tr <- 0
smol7$tr[smol7$three > -1] <- 1

smol8$tr <- 0
smol8$tr[smol8$three > -1] <- 1

smol9$tr <- 0
smol9$tr[smol9$four > -1] <- 1

smol10$tr <- 0
smol10$tr[smol10$four > -1] <- 1

smol11$tr <- 0
smol11$tr[smol11$five > -1] <- 1

smol12$tr <- 0
smol12$tr[smol12$five > -1] <- 1

smol13$tr <- 0
smol13$tr[smol13$six > -1] <- 1

smol14$tr <- 0
smol14$tr[smol14$six > -1] <- 1

smol15$tr <- 0
smol15$tr[smol15$seven > -1] <- 1

smol16$tr <- 0
smol16$tr[smol16$seven > -1] <- 1

smol17$tr <- 0
smol17$tr[smol17$eight > -1] <- 1

smol18$tr <- 0
smol18$tr[smol18$eight > -1] <- 1

smol19$tr <- 0
smol19$tr[smol19$nine > -1] <- 1

smol20$tr <- 0
smol20$tr[smol20$nine > -1] <- 1

smol21$tr <- 0
smol21$tr[smol21$ten > -1] <- 1

smol22$tr <- 0
smol22$tr[smol22$ten > -1] <- 1

smol23$tr <- 0
smol23$tr[smol23$eleven > -1] <- 1

smol24$tr <- 0
smol24$tr[smol24$eleven > -1] <- 1

smol25$tr <- 0
smol25$tr[smol25$twelve > -1] <- 1

smol26$tr <- 0
smol26$tr[smol26$twelve > -1] <- 1

smol27$tr <- 0
smol27$tr[smol27$thirteen > -1] <- 1

smol28$tr <- 0
smol28$tr[smol28$thirteen > -1] <- 1

smol29$tr <- 0
smol29$tr[smol29$thirteen > -1] <- 1

smol30$tr <- 0
smol30$tr[smol30$thirteen > -1] <- 1

smol31$tr <- 0
smol31$tr[smol31$fourteen > -1] <- 1

smol32$tr <- 0
smol32$tr[smol32$fourteen > -1] <- 1

smol33$tr <- 0
smol33$tr[smol33$fifteen > -1] <- 1

smol34$tr <- 0
smol34$tr[smol34$fifteen > -1] <- 1

smol35$tr <- 0
smol35$tr[smol35$fifteen > -1] <- 1

smol36$tr <- 0
smol36$tr[smol36$fifteen > -1] <- 1

smol37$tr <- 0
smol37$tr[smol37$run > -1] <- 1

smol38$tr <- 0
smol38$tr[smol38$run > -1] <- 1

smol39$tr <- 0
smol39$tr[smol39$run > -1] <- 1

smol40$tr <- 0
smol40$tr[smol40$run > -1] <- 1



# manual RDD

smol1$wt <- kernelwts(smol1$lamp, center = 0, bw = 30)
reg1 <- lm(factor ~ tr + lamp + tr*lamp, data = subset(smol1, wt > 0), weights = wt)
slope1 <- coeftest(reg1, vcov. = vcovHC(reg1, type = "HC1"))
coef1 <- slope1[4,1]
up1 <- slope1[4,1] + slope1[4,2]*1.96
low1 <- slope1[4,1] - slope1[4,2]*1.96

smol2$wt <- kernelwts(smol2$lamp, center = 0, bw = 30)
reg2 <- lm(factor ~ tr + lamp + tr*lamp, data = subset(smol2, wt > 0), weights = wt)
slope2 <- coeftest(reg2, vcov. = vcovHC(reg2, type = "HC1"))
## polynomials
reg2b <- lm(factor ~ tr + lamp + tr*lamp + I(lamp^2) + tr*I(lamp^2), data = subset(smol2, wt > 0), weights = wt)
slope2b <- coeftest(reg2b, vcov. = vcovHC(reg2b, type = "HC1"))
slope2b
##
coef2 <- slope2[4,1]
up2 <- slope2[4,1] + slope2[4,2]*1.96
low2 <- slope2[4,1] - slope2[4,2]*1.96

smol3$wt <- kernelwts(smol3$one, center = 0, bw = 30)
reg3 <- lm(factor ~ tr + one + tr*one + cntry, data = subset(smol3, wt > 0), weights = wt)
slope3 <- coeftest(reg3, vcov. = vcovHC(reg3, type = "HC1"))
coef3 <- slope3[13,1]
up3 <- slope3[13,1] + slope3[13,2]*1.96
low3 <- slope3[13,1] - slope3[13,2]*1.96

smol4$wt <- kernelwts(smol4$one, center = 0, bw = 30)
reg4 <- lm(factor ~ tr + one + tr*one + cntry, data = subset(smol4, wt > 0), weights = wt)
slope4 <- coeftest(reg4, vcov. = vcovHC(reg4, type = "HC1"))
## polynomials
reg4b <- lm(factor ~ tr + one + tr*one + I(one^2) + tr*I(one^2) + cntry, data = subset(smol4, wt > 0), weights = wt)
slope4b <- coeftest(reg4b, vcov. = vcovHC(reg4b, type = "HC1"))
slope4b
##
coef4 <- slope4[13,1]
up4 <- slope4[13,1] + slope4[13,2]*1.96
low4 <- slope4[13,1] - slope4[13,2]*1.96

smol5$wt <- kernelwts(smol5$two, center = 0, bw = 30)
reg5 <- lm(factor ~ tr + two + tr*two + cntry, data = subset(smol5, wt > 0), weights = wt)
slope5 <- coeftest(reg5, vcov. = vcovHC(reg5, type = "HC1"))
coef5 <- slope5[15,1]
up5 <- slope5[15,1] + slope5[15,2]*1.96
low5 <- slope5[15,1] - slope5[15,2]*1.96

smol6$wt <- kernelwts(smol6$two, center = 0, bw = 30)
reg6 <- lm(factor ~ tr + two + tr*two + cntry, data = subset(smol6, wt > 0), weights = wt)
slope6 <- coeftest(reg6, vcov. = vcovHC(reg6, type = "HC1"))
## polynomials
reg6b <- lm(factor ~ tr + two + tr*two + I(two^2) + tr*I(two^2) + cntry, data = subset(smol6, wt > 0), weights = wt)
slope6b <- coeftest(reg6b, vcov. = vcovHC(reg6b, type = "HC1"))
slope6b
##
coef6 <- slope6[15,1]
up6 <- slope6[15,1] + slope6[15,2]*1.96
low6 <- slope6[15,1] - slope6[15,2]*1.96

smol7$wt <- kernelwts(smol7$three, center = 0, bw = 30)
reg7 <- lm(factor ~ tr + three + tr*three + cntry, data = subset(smol7, wt > 0), weights = wt)
slope7 <- coeftest(reg7, vcov. = vcovHC(reg7, type = "HC1"))
coef7 <- slope7[18,1]
up7 <- slope7[18,1] + slope7[18,2]*1.96
low7 <- slope7[18,1] - slope7[18,2]*1.96

smol8$wt <- kernelwts(smol8$three, center = 0, bw = 30)
reg8 <- lm(factor ~ tr + three + tr*three + cntry, data = subset(smol8, wt > 0), weights = wt)
slope8 <- coeftest(reg8, vcov. = vcovHC(reg8, type = "HC1"))
## polynomials
reg8b <- lm(factor ~ tr + three + tr*three + I(three^2) + tr*I(three^2) + cntry, data = subset(smol8, wt > 0), weights = wt)
slope8b <- coeftest(reg8b, vcov. = vcovHC(reg8b, type = "HC1"))
slope8b
##
coef8 <- slope8[18,1]
up8 <- slope8[18,1] + slope8[18,2]*1.96
low8 <- slope8[18,1] - slope8[18,2]*1.96

smol9$wt <- kernelwts(smol9$four, center = 0, bw = 30)
reg9 <- lm(factor ~ tr + four + tr*four + cntry, data = subset(smol9, wt > 0), weights = wt)
slope9 <- coeftest(reg9, vcov. = vcovHC(reg9, type = "HC1"))
coef9 <- slope9[10,1]
up9 <- slope9[10,1] + slope9[10,2]*1.96
low9 <- slope9[10,1] - slope9[10,2]*1.96

smol10$wt <- kernelwts(smol10$four, center = 0, bw = 30)
reg10 <- lm(factor ~ tr + four + tr*four + cntry, data = subset(smol10, wt > 0), weights = wt)
slope10 <- coeftest(reg10, vcov. = vcovHC(reg10, type = "HC1"))
## polynomials
reg10b <- lm(factor ~ tr + four + tr*four + I(four^2) + tr*I(four^2) + cntry, data = subset(smol10, wt > 0), weights = wt)
slope10b <- coeftest(reg10b, vcov. = vcovHC(reg10b, type = "HC1"))
slope10b
##
coef10 <- slope10[10,1]
up10 <- slope10[10,1] + slope10[10,2]*1.96
low10 <- slope10[10,1] - slope10[10,2]*1.96

smol11$wt <- kernelwts(smol11$five, center = 0, bw = 30)
reg11 <- lm(factor ~ tr + five + tr*five + cntry, data = subset(smol11, wt > 0), weights = wt)
slope11 <- coeftest(reg11, vcov. = vcovHC(reg11, type = "HC1"))
coef11 <- slope11[15,1]
up11 <- slope11[15,1] + slope11[15,2]*1.96
low11 <- slope11[15,1] - slope11[15,2]*1.96

smol12$wt <- kernelwts(smol12$five, center = 0, bw = 30)
reg12 <- lm(factor ~ tr + five + tr*five + cntry, data = subset(smol12, wt > 0), weights = wt)
slope12 <- coeftest(reg12, vcov. = vcovHC(reg12, type = "HC1"))
## polynomials
reg12b <- lm(factor ~ tr + five + tr*five + I(five^2) + tr*I(five^2) + cntry, data = subset(smol12, wt > 0), weights = wt)
slope12b <- coeftest(reg12b, vcov. = vcovHC(reg12b, type = "HC1"))
slope12b
##
coef12 <- slope12[15,1]
up12 <- slope12[15,1] + slope12[15,2]*1.96
low12 <- slope12[15,1] - slope12[15,2]*1.96

smol13$wt <- kernelwts(smol13$six, center = 0, bw = 30)
reg13 <- lm(factor ~ tr + six + tr*six + cntry, data = subset(smol13, wt > 0), weights = wt)
slope13 <- coeftest(reg13, vcov. = vcovHC(reg13, type = "HC1"))
coef13 <- slope13[21,1]
up13 <- slope13[21,1] + slope13[21,2]*1.96
low13 <- slope13[21,1] - slope13[21,2]*1.96

smol14$wt <- kernelwts(smol14$six, center = 0, bw = 30)
reg14 <- lm(factor ~ tr + six + tr*six + cntry, data = subset(smol14, wt > 0), weights = wt)
slope14 <- coeftest(reg14, vcov. = vcovHC(reg14, type = "HC1"))
## polynomials
reg14b <- lm(factor ~ tr + six + tr*six + I(six^2) + tr*I(six^2) + cntry, data = subset(smol14, wt > 0), weights = wt)
slope14b <- coeftest(reg14b, vcov. = vcovHC(reg14b, type = "HC1"))
slope14b
##
coef14 <- slope14[21,1]
up14 <- slope14[21,1] + slope14[21,2]*1.96
low14 <- slope14[21,1] - slope14[21,2]*1.96

smol15$wt <- kernelwts(smol15$seven, center = 0, bw = 30)
reg15 <- lm(factor ~ tr + seven + tr*seven + cntry, data = subset(smol15, wt > 0), weights = wt)
slope15 <- coeftest(reg15, vcov. = vcovHC(reg15, type = "HC1"))
coef15 <- slope15[21,1]
up15 <- slope15[21,1] + slope15[21,2]*1.96
low15 <- slope15[21,1] - slope15[21,2]*1.96

smol16$wt <- kernelwts(smol16$seven, center = 0, bw = 30)
reg16 <- lm(factor ~ tr + seven + tr*seven + cntry, data = subset(smol16, wt > 0), weights = wt)
slope16 <- coeftest(reg16, vcov. = vcovHC(reg16, type = "HC1"))
## polynomials
reg16b <- lm(factor ~ tr + seven + tr*seven + I(seven^2) + tr*I(seven^2) + cntry, data = subset(smol16, wt > 0), weights = wt)
slope16b <- coeftest(reg16b, vcov. = vcovHC(reg16b, type = "HC1"))
slope16b
##
coef16 <- slope16[21,1]
up16 <- slope16[21,1] + slope16[21,2]*1.96
low16 <- slope16[21,1] - slope16[21,2]*1.96

smol17$wt <- kernelwts(smol17$eight, center = 0, bw = 30)
reg17 <- lm(factor ~ tr + eight + tr*eight + cntry, data = subset(smol17, wt > 0), weights = wt)
slope17 <- coeftest(reg17, vcov. = vcovHC(reg17, type = "HC1"))
coef17 <- slope17[21,1]
up17 <- slope17[21,1] + slope17[21,2]*1.96
low17 <- slope17[21,1] - slope17[21,2]*1.96

smol18$wt <- kernelwts(smol18$eight, center = 0, bw = 30)
reg18 <- lm(factor ~ tr + eight + tr*eight + cntry, data = subset(smol18, wt > 0), weights = wt)
slope18 <- coeftest(reg18, vcov. = vcovHC(reg18, type = "HC1"))
## polynomials
reg18b <- lm(factor ~ tr + eight + tr*eight + I(eight^2) + tr*I(eight^2) + cntry, data = subset(smol18, wt > 0), weights = wt)
slope18b <- coeftest(reg18b, vcov. = vcovHC(reg18b, type = "HC1"))
slope18b
##
coef18 <- slope18[21,1]
up18 <- slope18[21,1] + slope18[21,2]*1.96
low18 <- slope18[21,1] - slope18[21,2]*1.96

smol19$wt <- kernelwts(smol19$nine, center = 0, bw = 30)
reg19 <- lm(factor ~ tr + nine + tr*nine + cntry, data = subset(smol19, wt > 0), weights = wt)
slope19 <- coeftest(reg19, vcov. = vcovHC(reg19, type = "HC1"))
coef19 <- slope19[20,1]
up19 <- slope19[20,1] + slope19[20,2]*1.96
low19 <- slope19[20,1] - slope19[20,2]*1.96

smol20$wt <- kernelwts(smol20$nine, center = 0, bw = 30)
reg20 <- lm(factor ~ tr + nine + tr*nine + cntry, data = subset(smol20, wt > 0), weights = wt)
slope20 <- coeftest(reg20, vcov. = vcovHC(reg20, type = "HC1"))
## polynomials
reg20b <- lm(factor ~ tr + nine + tr*nine + I(nine^2) + tr*I(nine^2) + cntry, data = subset(smol20, wt > 0), weights = wt)
slope20b <- coeftest(reg20b, vcov. = vcovHC(reg20b, type = "HC1"))
slope20b
##
coef20 <- slope20[20,1]
up20 <- slope20[20,1] + slope20[20,2]*1.96
low20 <- slope20[20,1] - slope20[20,2]*1.96

smol21$wt <- kernelwts(smol21$ten, center = 0, bw = 30)
reg21 <- lm(factor ~ tr + ten + tr*ten + cntry, data = subset(smol21, wt > 0), weights = wt)
slope21 <- coeftest(reg21, vcov. = vcovHC(reg21, type = "HC1"))
coef21 <- slope21[15]
up21 <- slope21[15,1] + slope21[15,2]*1.96
low21 <- slope21[15,1] - slope21[15,2]*1.96

smol22$wt <- kernelwts(smol22$ten, center = 0, bw = 30)
reg22 <- lm(factor ~ tr + ten + tr*ten + cntry, data = subset(smol22, wt > 0), weights = wt)
slope22 <- coeftest(reg22, vcov. = vcovHC(reg22, type = "HC1"))
## polynomials
reg22b <- lm(factor ~ tr + ten + tr*ten + I(ten^2) + tr*I(ten^2) + cntry, data = subset(smol22, wt > 0), weights = wt)
slope22b <- coeftest(reg22b, vcov. = vcovHC(reg22b, type = "HC1"))
slope22b
##
coef22 <- slope22[15]
up22 <- slope22[15,1] + slope22[15,2]*1.96
low22 <- slope22[15,1] - slope22[15,2]*1.96

smol23$wt <- kernelwts(smol23$eleven, center = 0, bw = 30)
reg23 <- lm(factor ~ tr + eleven + tr*eleven + cntry, data = subset(smol23, wt > 0), weights = wt)
slope23 <- coeftest(reg23, vcov. = vcovHC(reg23, type = "HC1"))
coef23 <- slope23[8]
up23 <- slope23[8,1] + slope23[8,2]*1.96
low23 <- slope23[8,1] - slope23[8,2]*1.96

smol24$wt <- kernelwts(smol24$eleven, center = 0, bw = 30)
reg24 <- lm(factor ~ tr + eleven + tr*eleven + cntry, data = subset(smol24, wt > 0), weights = wt)
slope24 <- coeftest(reg24, vcov. = vcovHC(reg24, type = "HC1"))
## polynomials
reg24b <- lm(factor ~ tr + eleven + tr*eleven + I(eleven^2) + tr*I(eleven^2) + cntry, data = subset(smol24, wt > 0), weights = wt)
slope24b <- coeftest(reg24b, vcov. = vcovHC(reg24b, type = "HC1"))
slope24b
##
coef24 <- slope24[8]
up24 <- slope24[8,1] + slope24[8,2]*1.96
low24 <- slope24[8,1] - slope24[8,2]*1.96

smol25$wt <- kernelwts(smol25$twelve, center = 0, bw = 30)
reg25 <- lm(factor ~ tr + twelve + tr*twelve + cntry, data = subset(smol25, wt > 0), weights = wt)
slope25 <- coeftest(reg25, vcov. = vcovHC(reg25, type = "HC1"))
coef25 <- slope25[7]
up25 <- slope25[7,1] + slope25[7,2]*1.96
low25 <- slope25[7,1] - slope25[7,2]*1.96

smol26$wt <- kernelwts(smol26$twelve, center = 0, bw = 30)
reg26 <- lm(factor ~ tr + twelve + tr*twelve + cntry, data = subset(smol26, wt > 0), weights = wt)
slope26 <- coeftest(reg26, vcov. = vcovHC(reg26, type = "HC1"))
## polynomials
reg26b <- lm(factor ~ tr + twelve + tr*twelve + I(twelve^2) + tr*I(twelve^2) + cntry, data = subset(smol26, wt > 0), weights = wt)
slope26b <- coeftest(reg26b, vcov. = vcovHC(reg26b, type = "HC1"))
slope26b
##
coef26 <- slope26[7]
up26 <- slope26[7,1] + slope26[7,2]*1.96
low26 <- slope26[7,1] - slope26[7,2]*1.96

smol27$wt <- kernelwts(smol27$thirteen, center = 0, bw = 30)
reg27 <- lm(factor ~ tr + thirteen + tr*thirteen + cntry, data = subset(smol27, wt > 0), weights = wt)
slope27 <- coeftest(reg27, vcov. = vcovHC(reg27, type = "HC1"))
coef27 <- slope27[24]
up27 <- slope27[24,1] + slope27[24,2]*1.96
low27 <- slope27[24,1] - slope27[24,2]*1.96

smol28$wt <- kernelwts(smol28$thirteen, center = 0, bw = 30)
reg28 <- lm(factor ~ tr + thirteen + tr*thirteen + cntry, data = subset(smol28, wt > 0), weights = wt)
slope28 <- coeftest(reg28, vcov. = vcovHC(reg28, type = "HC1"))
## polynomials
reg28b <- lm(factor ~ tr + thirteen + tr*thirteen + I(thirteen^2) + tr*I(thirteen^2) + cntry, data = subset(smol28, wt > 0), weights = wt)
slope28b <- coeftest(reg28b, vcov. = vcovHC(reg28b, type = "HC1"))
slope28b
##
coef28 <- slope28[23,1]
up28 <- slope28[23,1] + slope28[23,2]*1.96
low28 <- slope28[23,1] - slope28[23,2]*1.96

smol29$wt <- kernelwts(smol29$thirteen, center = 0, bw = 30)
reg29 <- lm(factor ~ tr + thirteen + tr*thirteen, data = subset(smol29, wt > 0), weights = wt)
slope29 <- coeftest(reg29, vcov. = vcovHC(reg29, type = "HC1"))
coef29 <- slope29[4,1]
up29 <- slope29[4,1] + slope29[4,2]*1.96
low29 <- slope29[4,1] - slope29[4,2]*1.96

smol30$wt <- kernelwts(smol30$thirteen, center = 0, bw = 30)
reg30 <- lm(factor ~ tr + thirteen + tr*thirteen, data = subset(smol30, wt > 0), weights = wt)
slope30 <- coeftest(reg30, vcov. = vcovHC(reg30, type = "HC1"))
## polynomials
reg30b <- lm(factor ~ tr + thirteen + tr*thirteen + I(thirteen^2) + tr*I(thirteen^2), data = subset(smol30, wt > 0), weights = wt)
slope30b <- coeftest(reg30b, vcov. = vcovHC(reg30b, type = "HC1"))
slope30b
##
coef30 <- slope30[4,1]
up30 <- slope30[4,1] + slope30[4,2]*1.96
low30 <- slope30[4,1] - slope30[4,2]*1.96

smol31$wt <- kernelwts(smol31$fourteen, center = 0, bw = 30)
reg31 <- lm(factor ~ tr + fourteen + tr*fourteen + cntry, data = subset(smol31, wt > 0), weights = wt)
slope31 <- coeftest(reg31, vcov. = vcovHC(reg31, type = "HC1"))
coef31 <- slope31[10,1]
up31 <- slope31[10,1] + slope31[10,2]*1.96
low31 <- slope31[10,1] - slope31[10,2]*1.96

smol32$wt <- kernelwts(smol32$fourteen, center = 0, bw = 30)
reg32 <- lm(factor ~ tr + fourteen + tr*fourteen + cntry, data = subset(smol32, wt > 0), weights = wt)
slope32 <- coeftest(reg32, vcov. = vcovHC(reg32, type = "HC1"))
## polynomials
reg32b <- lm(factor ~ tr + fourteen + tr*fourteen + I(fourteen^2) + tr*I(fourteen^2) + cntry, data = subset(smol32, wt > 0), weights = wt)
slope32b <- coeftest(reg32b, vcov. = vcovHC(reg32b, type = "HC1"))
slope32b
##
coef32 <- slope32[10,1]
up32 <- slope32[10,1] + slope32[10,2]*1.96
low32 <- slope32[10,1] - slope32[10,2]*1.96

smol33$wt <- kernelwts(smol33$fifteen, center = 0, bw = 30)
reg33 <- lm(factor ~ tr + fifteen + tr*fifteen + cntry, data = subset(smol33, wt > 0), weights = wt)
slope33 <- coeftest(reg33, vcov. = vcovHC(reg33, type = "HC1"))
coef33 <- slope33[16,1]
up33 <- slope33[16,1] + slope33[16,2]*1.96
low33 <- slope33[16,1] - slope33[16,2]*1.96

smol34$wt <- kernelwts(smol34$fifteen, center = 0, bw = 30)
reg34 <- lm(factor ~ tr + fifteen + tr*fifteen + cntry, data = subset(smol34, wt > 0), weights = wt)
slope34 <- coeftest(reg34, vcov. = vcovHC(reg34, type = "HC1"))
## polynomials
reg34b <- lm(factor ~ tr + fifteen + tr*fifteen + I(fifteen^2) + tr*I(fifteen^2) + cntry, data = subset(smol34, wt > 0), weights = wt)
slope34b <- coeftest(reg34b, vcov. = vcovHC(reg34b, type = "HC1"))
slope34b
##
coef34 <- slope34[16,1]
up34 <- slope34[16,1] + slope34[16,2]*1.96
low34 <- slope34[16,1] - slope34[16,2]*1.96

smol35$wt <- kernelwts(smol35$fifteen, center = 0, bw = 30)
reg35 <- lm(factor ~ tr + fifteen + tr*fifteen, data = subset(smol35, wt > 0), weights = wt)
slope35 <- coeftest(reg35, vcov. = vcovHC(reg35, type = "HC1"))
coef35 <- slope35[4,1]
up35 <- slope35[4,1] + slope35[4,2]*1.96
low35 <- slope35[4,1] - slope35[4,2]*1.96

smol36$wt <- kernelwts(smol36$fifteen, center = 0, bw = 30)
reg36 <- lm(factor ~ tr + fifteen + tr*fifteen, data = subset(smol36, wt > 0), weights = wt)
slope36 <- coeftest(reg36, vcov. = vcovHC(reg36, type = "HC1"))
## polynomials
reg36b <- lm(factor ~ tr + fifteen + tr*fifteen + I(fifteen^2) + tr*I(fifteen^2), data = subset(smol36, wt > 0), weights = wt)
slope36b <- coeftest(reg36b, vcov. = vcovHC(reg36b, type = "HC1"))
slope36b
##
coef36 <- slope36[4,1]
up36 <- slope36[4,1] + slope36[4,2]*1.96
low36 <- slope36[4,1] - slope36[4,2]*1.96

smol37$wt <- kernelwts(smol37$run, center = 0, bw = 30)
reg37 <- lm(factor ~ tr + run + tr*run + cntry, data = subset(smol37, wt > 0), weights = wt)
slope37 <- coeftest(reg37, vcov. = vcovHC(reg37, type = "HC1"))
coef37 <- slope37[35,1]
up37 <- slope37[35,1] + slope37[35,2]*1.96
low37 <- slope37[35,1] - slope37[35,2]*1.96

smol38$wt <- kernelwts(smol38$run, center = 0, bw = 30)
reg38 <- lm(factor ~ tr + run + tr*run, data = subset(smol38, wt > 0), weights = wt)
slope38 <- coeftest(reg38, vcov. = vcovHC(reg38, type = "HC1"))
coef38 <- slope38[4,1]
up38 <- slope38[4,1] + slope38[4,2]*1.96
low38 <- slope38[4,1] - slope38[4,2]*1.96

smol39$wt <- kernelwts(smol39$run, center = 0, bw = 30)
reg39 <- lm(factor ~ tr + run + tr*run + cntry, data = subset(smol39, wt > 0), weights = wt)
slope39 <- coeftest(reg39, vcov. = vcovHC(reg39, type = "HC1"))
coef39 <- slope39[35,1]
up39 <- slope39[35,1] + slope39[35,2]*1.96
low39 <- slope39[35,1] - slope39[35,2]*1.96

smol40$wt <- kernelwts(smol40$run, center = 0, bw = 30)
reg40 <- lm(factor ~ tr + run + tr*run, data = subset(smol40, wt > 0), weights = wt)
slope40 <- coeftest(reg40, vcov. = vcovHC(reg40, type = "HC1"))
coef40 <- slope40[4,1]
up40 <- slope40[4,1] + slope40[4,2]*1.96
low40 <- slope40[4,1] - slope40[4,2]*1.96

summary(reg38)

# plot slopes and CI

load("multi.RData")

multi$coef <- c(coef1, coef3,coef5,coef7,coef9,coef11,coef13,coef15,coef17,
                coef19,coef21,coef23,coef25,coef27,coef29, coef31, coef33, 
                coef35,coef37,coef38)

multi$lowerci <- c(low1, low3,low5,low7,low9,low11,low13,low15,low17,low19,
                   low21,low23,low25,low27, low29,low31,low33,low35,
                   low37,low38)

multi$upperci <- c(up1, up3,up5,up7,up9,up11,up13,up15,up17,up19,
                   up21,up23,up25,up27, up29, up31, up33, up35,
                   up37,up38)



multi$date <- factor(multi$date, levels=unique(multi$date))
appendixs <- ggplot(multi, aes(x = date, y = coef, color = Sample)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 7, position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin = lowerci, ymax = upperci), size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-0.5,0.5) +
  scale_colour_manual(values = c("Italy only" = "#F8766D", "Cross-national" = "#00BFC4")) +
  theme_light(base_size = 30) +
  scale_x_discrete(guide = guide_axis(angle = 90)) +
  xlab("Date of the Event") +
  ylab("RDD Slope Change")


multi$coefb <- c(coef2, coef4,coef6,coef8,coef10,coef12,coef14,coef16,coef18,
                 coef20,coef22,coef24,coef26,coef28, coef30, coef32, coef34, coef36,
                 coef39,coef40)

multi$lowercib <- c(low2, low4,low6,low8,low10,low12,low14,low16,low18,low20,
                    low22,low24,low26,low28, low30, low32, low34, low36,
                    low39, low40)

multi$uppercib <- c(up2, up4,up6,up8,up10,up12,up14,up16,up18,up20,
                    up22,up24,up26,up28,up30, up32, up34, up36,
                    up39,up40)


multi$date <- factor(multi$date, levels=unique(multi$date))
mains <- ggplot(multi, aes(x = date, y = coefb, color = Sample)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed", size = 1) +
  geom_point(size = 7, position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin = lowercib, ymax = uppercib), size = 2, width = 0.2, position = position_dodge(0.5)) +
  ylim(-0.2,0.2) +
  scale_colour_manual(values = c("Italy only" = "#F8766D", "Cross-national" = "#00BFC4")) +
  theme_light(base_size = 30) +
  scale_x_discrete(guide = guide_axis(angle = 90)) +
  xlab("Date of the Event") +
  ylab("RDD Slope Change")

ggarrange(mainc, mains)
ggarrange(appendixc, appendixs)
